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We present a formalism for computing the complexity of metastable states and the zero- 
temperature magnetic hysteresis loop in the soft-spin random-field model in finite dimensions. The 
complexity is obtained as the Legendre transform of the free-energy associated to a certain ac- 
tion in replica space and the hysteresis loop above the critical disorder is defined as the curve in the 
field-magnetization plane where the complexity vanishes; the nonequilibrium magnetization is there- 
fore obtained without having to follow the dynamical evolution. We use approximations borrowed 
from condensed-matter theory and based on assumptions on the structure of the direct correlation 
functions (or proper vertices), such as a local approximation for the self-energies, to calculate the 
hysteresis loop in three dimensions, the correlation functions along the loop, and the second moment 
of the avalanche-size distribution. 
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I. INTRODUCTION 

The random-field Ising model (RFIM) at zero temperature is a prototype of many disordered systems which exhibit 
hysteretic and jerky behavior when an external parameter (e.g. magnetic field, pressure, strain) is slowly changed[T|. 
This behavior is related to the existence of a corrugated (free) energy landscape with many local minima (or metastable 
states) separated by barriers much larger than fc^T, therefore preventing relaxation towards equilibrium on exper- 
imental time scales. As a consequence, the response to an external driving field is a series of discontinuous jumps 
(called avalanches) from one metastable state to another, the number and size of these jumps varying with the amount 
of disorder. In the 3-dimensional RFIM with Gaussian random fields two regimes of avalanches are observed [2 [3]: at 
large disorder, there are many microscopic jumps resulting in a smooth magnetization curve macroscopically; at low 
disorder, there is a system-spanning avalanche resulting in a macroscopic jump in the magnetization curve. These 
two regimes are separated by a critical point at which avalanches of all sizes are observed. As shown recently, this 
type of nonequilibrium disorder-induced phase transition may underly the hysteretic behavior of 4 He adsorbed in high 
porosity silica aerogels [H |5]. 

Even though the T = RFIM and similar models have been extensively studied, there is currently no analytical 
tool to compute the saturation hysteresis loop in finite dimensions (even approximately) and no theory to describe 
the statistical properties of the metastable states, for instance their configurational entropy (also called 'complexity') 
as a function of magnetization. Furthermore, the behavior of the correlation functions along the hysteresis loop is 
unknown, although this is an issue of experimental interest [6]. Of course, describing analytically such a nonequilibrium 
situation where the present state of the system depends on its past history is a difficult problem and at first sight 
there seems to be no other choice than to follow the dynamical evolution for given initial conditions (e.g. one of the 
two saturated states corresponding to infinitely positive or negative magnetic field) . In this way, one can treat exactly 
'mean-field' systems, i.e. fully connected lattices [2J |3] or lattices with a locally tree-like structure such as the Bethc 
lattice [5], and obtain analytical expressions for the saturation hysteresis loop or the avalanche size distribution. To 
go further and build a field-theoretical description of these phenomena it then seems necessary to employ the Martin- 
Siggia-Rose formalism [3 . In ferromagnetic systems, however, an alternative route is possible thanks to a remarkable 
property of the saturation hysteresis loop in the large-disorder regime: the loop is the convex hull of the metastable 
states in the two-dimensional field-magnetization plane |5Hl2|. In other words, the complexity is zero outside the loop 
and positive inside [13]. Moreover, it exactly vanishes along the loop since there is a single 'extremal' metastable state 
at a given field (at least for a continuous distribution of the random field). Determining the hysteresis loop is then 
tantamount to counting the number of metastable states at fixed field as a function of their magnetization and finding 
the magnetization at which the complexity vanishes. The main difficulty is that one must count the typical (i.e. most 
likely) number of metastable states and compute the associated quenched complexity, which of course is a nontrivial 
task requiring the use of replicas. Nevertheless, at least in principle, one can thereby study hysteresis and avalanche 
statistics without referring to the dynamical evolution. 

In the present paper we test this approach by studying the soft-spin version of the RFIM introduced in Ref.[3], 
where the spins take on continuous values between — oo and +oo and are confined by a bistable potential with a 
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linear cusp. This model has the advantage over the standard RFIM to show hysteresis for all values of the disorder 
in mean-field theory and is anyhow expected to behave as the hard-spin model on long length scales. Our basic 
strategy is to express the complexity as the Legendre transform of the 'free-energy' associated to a certain action in 
replica space and then to use approximations borrowed from condensed-matter theory to calculate the corresponding 
correlation (or Green's) functions and thermodynamic properties. Our focus on the correlation functions is also 
motivated by analytical and numerical calculations on the Bethe lattice which suggest that their spatial structure 
along the hysteresis loop is simpler than at equilibrium|15 . 

The paper is organized as follows. In section II, we define the model and the observables that we want to compute. In 
section III, we introduce the general formalism in the context of a replica method, focusing on the various correlation 
functions. In section IV, we first consider the random-phase approximation (RPA) which is equivalent to mean- 
field theory and becomes exact in infinite dimension. In section V, we go beyond the RPA by introducing a local 
self-energy approximation (LSEA) and we compare the predictions for the magnetization, the correlation functions, 
and the second-moment of the avalanche-size distribution along the hysteresis loop to simulation data. Concluding 
remarks and directions for future work are provided in Section VI. Additional details on the analytical calculations 
are provided in the appendix. 



II. MODEL AND OBSERVABLES 



We consider a collection of N soft spins placed on the sites of a d-dimensional hypercubic lattice and interacting 
via the Hamitonian 

h = - j s ^ - T,( H + h *>* + E v ^ w 

<i,j> i i 

where J > is a ferromagnetic interaction that couples nearest-neighbor spins, H is an external uniform field, and 
{hi} is a set of quenched random fields drawn independently from a Gaussian probability distribution p(h) with zero 
mean and variance A. V(s) is a double- well potential for which we choose the same cuspy form as in Ref . pfl IT2"] 

y(s)^[ S -sign( S )] 2 (2) 

so that all solutions of dH/dsi = 0, i.e. all stationary points of the Hamiltonian, are local minima. This greatly 
facilitates the present study as there is no need to explicitly discard local maxima and saddle-points through the 
consideration of the Hessian of the energy function. These local minima are the so-called 'metastable' states in which 
each spin satisfies 



H + JE 



s l - sign(s l ) = ^ 3/l : (3) 

k 

where j is a neighbor of i on the lattice. 

At zero-temperature, one can define a local relaxation dynamics in which each spin is forced to satisfy Eq. [3] as 
the external field is changed[3]. When adiabatically varying H from — oo to +oo and backwards, the model exhibits 
hysteresis and the magnetization typically behaves as shown in Fig. 1. (The figure displays the results of a single 
simulation on a cubic lattice with k — 8 and A = 4|16j.) These values are arbitrarily chosen and will stay fixed in 
the rest of this work. The shape of the loop then changes with the coupling J: one can see that J = 0.3 and J = 0.6 
correspond to the large- and small-disorder regimes respectively, with a macroscopic jump in the latter case. 

For a given realization of the disorder, i.e. a set of random fields h = {hi}, and a given value of the external uniform 
field H, each metastable state is characterized at a macroscopic level by its magnetization, its energy, etc. In the 
present study, we only focus on the magnetization since this is sufficient to unambiguously determine the hysteresis 
loop. Our goal is then to compute properties averaged over all metastable states with a given magnetization m per site 
at a given external magnetic field H (with a flat measure) . A central quantity is the quenched complexity Eg (to, H ) , 
which is defined as 

f 
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E Q (m,#)= lim -lnAf(m,#;h) (4) 

N— >co iv 

where TV" (to, H; h) is the number of metastable states with magnetization to at the field H and the overbar denotes 
an average over the random-field distribution. We are also interested in the following two-point correlation (Green's) 
functions: 
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FIG. 1: Magnetization curves for the nonequilibrium zero-temperature random-field soft-spin model on a cubic lattice for 
k — 8, A = 4, J — 0.3 (black) and J = 0.6 (red). The simulation data correspond to a single disorder realization of linear size 
L = 100 using an increment in the external field SH = 0.01. 



(i) the spin-spin correlation functions, 

G SStij {m,H) = Gij(m,H) = G Ctij (m,H) + G d ,ij(m,H), (5) 
with the (disorder) connected and disconnected components defined as 

G c>i j (m, H) = < s t Sj > - < Si >< sj > (6) 

and 

Gd,ij{m, H) = < Si >< Sj > - < Si > < Sj > , (7) 

the brackets denoting an averaged over all metastable states with magnetization m at the field H , 

(ii) the correlation function involving the spin variable and the random field, 

G sh ,ij{m,H) = < Si> hj , (8) 

(iii) the spin-spin correlation function for two copies of the same disordered system coupled to different external 
fields and with different magnetizations, 

G d ,ij(m a , H a ; m b , H b ) = < Si > a < s 3 > b - < Si > a < Sj > b , (9) 

where the subscrit a indicates that {m a , H a } are fixed and similarly for the subscript b. 

Along the hysteresis loop, the complexity is zero (at least in the large disorder regime where the loop is continuous) [SJ- 
I12j and the connected spin-spin correlation is identically zero, as the fluctuations inside a metastable state vanish at 
zero temperature. On the other hand, the disconnected spin-spin correlation function and the spin-random-field one 
should be nontrivial. 

Finally, we are also interested in the distribution of avalanche sizes along the loop (say, along the ascending branch) . 
For a given disorder realization, the magnetization curve m(H;h) — (1/iV) consists of smooth parts resulting 

from the harmonic response to the magnetic field and a series of jumps of size S a (h) occurring at the fields H a (h) 
(these jumps are not visible in Fig. 1 due to the scale of the figure) [15] . The magnetization can thus be decomposed 
as 

m(H; h) = m smooth (H; h) + ^ SJ{H - H a ) (10) 

OL 

where 9(x) is the Heaviside step function. From this, one defines the jump (avalanche) density 



p(S,H) = J2 S (S - S a )S(H - H a ), 

a 



(ii) 
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so that p(S, H)dSdH is the number of avalanches of size between S and S + dS when the field is increased from H to 
H + dH (note that S is measured per site). In present work we will only compute the 'unnormalized' second moment 
/ S 2 p(S,H)dS — J2 a Sa5(H — H a ) (it is 'unnormalized' because p(S,H) as such is not a probability density and 
J p(S, H)dS is the total number of avalanches between H and H + dH). 



III. FORMALISM 

In order to control the local magnetization, we introduce an additional source H = {Hi} that is linearly coupled 
to the spins s = {si} and we consider the following (disorder-dependent) 'partition function' in an external magnetic 
field which is momentarily taken as nonuniform, H = {Hi}: 

ZI H,H;h] = /ps^. n <f) 

i 

= [vse^HSiVis^-jJ^Sj-Hi-hi] (12) 

1 j/i 

where the symbol Us refers to the integration over all the spin variables, Us = dsi...dsfj, and no Jacobian is needed 
(it would merely introduce a constant multiplicative factor). This partition function can be put into a more standard 
form by replacing the Dirac delta function by its Fourier representation, 

Z[H, H; h] = J UsUs e -S[M,h]+ft. B +H.a ( 13 ) 

where the action S is defined by 

S[s,s,h] =^2s i \V'{s i )- JY, S J-H (14) 

i j ft 

and s = {£i} are auxiliary (imaginary) variables; for conciseness, the factor l/(2i7r) associated to the integration of 
§i along the imaginary axis is adsorbed into d§i. 

This defines a 'free energy' VF[H,H;h] = lnZ[H, H;h] which is a random object whose cumulants give access to 
full information about the system, including the complexity and the correlation functions. The complexity is obtained 
from the first cumulant, 



W 1 [H,H]=W[U,H;h}, (15) 

via a Legendre transform, where 

ra ,[H,H] = a «, (16) 
a Hi 

and which for uniform sources takes the form 

X Q {m,H) = ±W 1 (H,H)-m(H,H)H. (17) 

(Here and below we use square brackets [...] when the arguments are locally varying and parenthesis (...) when they 
are uniform.) As discussed in Refs. 9-12 , the hysteresis loop in the large-disorder regime identifies with the curve 
£q(?ti, H) = in the limit H — > ±oo whereas the typical properties of the metastable states are obtained for H = 
which corresponds to the maximum of the complexity. 

The information about the distribution of avalanche sizes is contained in the higher-order cumulants 
W2P 1 , H 1 ; H 2 , H 2 ], iy 3 [H\H^H 2 ,H 2 ;H 3 ,H 3 ],..., where 



Vl^H 1 , H 1 ; H 2 , H 2 ] = W[H\ H 1 ; h]lF[H 2 , H 2 ; h] - W^H 1 , H 1 ; h] IF[H 2 ,H 2 ;h] (18) 



etc. 
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The correlation (Green's) functions are obtained by derivation of the cumulants with respect to the sources. For 
instance, the physical two-point spin-spin correlation functions introduced above are given by 

G ct3 = dW ^\ (19) 
c,lJ OHidHj 



dVt^tH 1 , H 1 , H 2 , H 2 ] 
- OH 1 dH? Ihi=h2=h,hi=h2=h » \ M ) 

where both right-hand sides are evaluated for uniform sources and H is considered as a function of m and H through 
the Legendre transform in Eq. ( |17[ ). The spin-random-field correlation function requires a little more thought. By 
using the property of Gaussian distributions, 



J dhp(h)hA(h) = -&J dh^^-A{h) = A J dh P( h ) 



dA{h) 
dh ' 



(21) 



one finds 



<s^h j = d ^m hj = A = A G cij , (22) 

J dHi J dHidHj * 3 V ; 



where 



Gc ^ ~ dHtdHj • (23) 
When all the sources are uniform, this yields 

dm(H,H) » . , 11 — 

dH = G ^ = °) = a N £ < Si > hj ■ (24) 

where G c (k) is the Fourier transform of G c ^. (In particular, this equation is valid in the limit H — > ±oo, i.e. along 
the hysteresis loop: surprisingly, it seems that this extension of the 'susceptibily sum-rule' to the nonequilibrium 
magnetization curve has not been noticed before.) 

To compute the average over disorder, the common procedure is to replicate the system n times and take the limit 
n — > at the end. (This is in contrast with the Martin-Siggia-Rose formalism in which the partition function is directly 
averaged over disorder, making the use of replicas unnecessary [3]; here indeed, the partition function in Eq. (13 1 is 
nontrivial so that one must average hiZ and not simply Z.) If one is interested in computing the cumulants of the 
random free-energy for generic arguments, one must introduce sources that act separately on each replica [T9U21j . 
After performing the average over the random- field distribution, we obtain a 'replica partition function', 

~ n 

Z rep [{H a ,H Q }] = / \[Vs a Vs a e -S™p[{sV°}]+£JH<\s°+H<\§°] ( 25 ) 

^ a=l 

with the replicated action given by 

s rep [{s\ s»}] = E E s ti v '( s t) J E s j] - j E E § " § i ■ ( 26 ) 

i a j/i i a,b 

The 'thermodynamic potential' W re p[{H a , H a }] = In Z rep [{H a , H a }] can then be expanded in increasing number of 
free replica sums [T5H2T] , 

oo 1 

W rep [{U a ,H a }} =E"| E W / p [H ai ,H ai ;H a2 ,H a2 ;...;H^,H a "] (27) 

p=l " a 1 ,a 2 ,..,a p 



G 



where the Wp's are continuous and symmetric functions of their arguments. This coincides with the cumulant 
expansion [20 j. The thermodynamic potential W rep [{H a , H a }] generates the 'magnetizations' mf and mf 

^ =< sf >= mf 

BE* 
BW 

dW ^-<sf>=mf, (28) 



and the correlation (or Green's) functions, e.g. at the pair level, 

G$ = dWre l 1 =< sfs^ >-<sf>< s* > 
3 BmBH" 3 3 

1 J 

6$ = dWre l =< sft >-<sl >< § b > 

13 dHi, a dHj, b 13 13 

£,ab _ 9W rep |„|6 _ Ja ><c |6 > (on) 

ij ~ 3H a BH b 1 J ' ^ ' 

where < ... > denotes an average over the replicated action. As ^^[{H", H a }] above, the magnetizations and the 
correlation functions can be expanded in increasing number of free replica sums: 

m l!a [{H e , H e }] = mf 1 [H a , H a ] + ]T mf 1 [H a , H a |H e , H e ] + \ £ mf 1 [H a , H a |H e , H e ; H^, ft?] + (30) 

and similarly for mf, whereas after decomposing the two-point functions as Gf b = Gf^Sab + ( wnere ^aHj ^ oes 
not contain any explicit Kronecker delta) , one has [20j [21] 



G^-[{H e ,H e }] = G^\j [H a , H a ] + ^c\j P", H a |H e , H e ] ^ £ G^- H a |^ H e ; H', + (31) 



2 

e,/ 



G^[{H e ,H e }] = G[%[^H a ;H fc ,H b ] + ^G^.[H a ,H a ;H fc ,H b |H e ,H e ] 

e 

+ \ E G Si [H a , H a ;H h ,H b |H e ,H e ;H / ,H- f ] + (32) 

e,/ 

and similarly for G«, 3 , G^ and 6^, G% y 

In the present approach, however, the central object is not W rep but the 'effective action' T rep which is the Legendrc 
transform of W rep with respect to the two sets of sources {H a } and {H a }, 

r rep [{m a , rh a }] = -W rep [{H a , H a }] + £ (m a .H Q + m Q .H a ) , (33) 

a 

so that 

BT 

jja _ UL rep 

1 ~ Bmf 
BT 

= (34) 
Bmf y ' 

T rep is the generating functional of the 'direct' correlation functions [or one-particle irreducible (1PI) functions, or 
else proper vertices, in field-theoretic language]. At the pair level, 

(jab _ 9^rep 



lJ BmfBmh 

(Jab _ 9T rep 

y BmfBrh b 



j 

(jab _ 9T rep 
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As above, T rep , H?,Hf, Cff, etc... can be expanded in increasing number of free replica sums. 

In the following, all the two-point functions will be put together as components of a matrix with both replica indices 
and spatial coordinates 



G = 




(36) 



where G, G, G have for elements Gf b , Gf b , Gf b ; a similar notation is used for the direct correlation functions, all 
collected in C. The matrix C is then just the inverse of G, i.e. 

q = or 1 . (37) 

The matrix inversion with respect to spatial coordinates is easily realized by Fourier transformation when all the 
sources are taken as uniform. The inversion with respect to the replica indices can be performed by using the expansion 
in free replica sums for G and C (see Eqs. ( 31|32 )) and proceeding to a term-by-term identification. The zeroth-order 
terms are then given by 



£M(k;m a ,m a ) = d 0] (k; m a , to")" 1 

G = f(k;m a ,m a ;m b 7 m b ) = -d 0] (k; m a , m'j-^fk; m a , m a ; m h , m b )Q l ° ] (k; m b , m 6 )" 1 , 



(38) 
(39) 



where Gj ^ are 2x2 matrices containing the components G Cj£ j, G c , dl G c ^ d as in Eq. (pL and similarly for C^. Note 



that the zeroth-order functions obtained when fixing the replica sources {H a ,H a } and fixing the replica 'magneti- 
zations' {m a ,m a } are generically related through G(H a , H a ) = G{m^(H a ,H a ),m^(H a ,H a )) with m' ' and to™ 
being the zeroth-order expressions as in Eq. (28). Since the zeroth-order contributions already contain the physics of 
the problem, we will not consider higher-order terms and will drop the superscript [0] in the following. 

When all replica sources or magnetizations are taken as equal, and provided that this limiting process is regular 
enough (this will be discussed in more detail below), replica symmetry is recovered and one obtains a set of 'Ornstein- 
Zernike' (OZ) equations (in the language of liquid-state theory [22]) which takes the following explicit form: 



G c (k) = 



G c (k) 



G c (k) = - 



G c (k)G c (k) - G c (k) 2 
G c (k) 



G c (k) = 



G c (k)G c (k) - G c (k) 2 
C c (k) 

c c (k)4(k) - c c (k) 2 



(40) 



and 



G d (k) = -[C d (k)G c (k) 2 + 2(7 d (k)G c (k)G c (k) + C d (k)G c (k) 2 } 

G d (k) = -[C d (k)G c (k)G c (k) + C d (k)(G c (k)6 c (k) + G c (k) 2 ) + 6 d (k)6 c (k)C c (k)} 

6 d (k) = -[d d (k)d c (kj 2 + 2G d (k)4(k)G c (k) + G d (k)G c (k) 2 ] , 



(41) 



where all functions depend on m and m. 

When the replica sources are not equal, the zeroth-order terms have the capability to capture a nonanalytic behavior 
in the dependence on the replica sources or replica magnetizations. This important feature is already displayed by 
the noninteracting system cor respo nding to J = (hereafter called the 'reference' system) whose properties are 

listed in the appendix (Eq. (A13) and below with K r — and K d = A). For instance, the 2-replica function 
G d ,ij(H a , H; H b , H) behaves like 



G d ,ij(H a , H; H b , H) = G d ^j(H, H; H, H) 



G d yf(H,H)\H a 



H"\ + 0((H a - H b ) z ) 



(42) 



when H a , H b — > H . We stress that this is not a spurious behavior due to the presence of a cusp in V(s). Indeed, even 
if V(s) were a smooth double- well potential, it would be necessary to distinguish between the two minima in order to 
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count the metastable states and, whatever the method, this would introduce a nonanalyticity in VF[H, H; h] and lead 
to a cusp in the 2-replica function Gd- This behavior is thus an intrinsic feature of the problem under consideration 
and is also intimately connected to the magnetization discontinuities along the hysteresis loop, i.e. to the avalanches. 
A similar connection is discussed in Rcf. [23 in the context of random elastic systems where the statistics of static 
avalanches (or 'shocks') is studied via the functional renormalization group. (In this case, however, the cusp only 
shows up in the course of the renormalization flow whereas it is always present here, even in the large disorder regime.) 

To show that the cusp is related to the (unnormalized) second moment of the avalanche size distribution, we consider 
the quantity 

G d (k = 0; H a ;H b ) = N \m(H a ; h)m(H b ; h) - m(H a )m(H b )] , (43) 
where it is implicit that the lim it H — > ±oo corresponding to the hysteresis loop (in the large-disorder regime) has 



been taken. Then, using Eqs. (10 1 and (42 1, taking the second derivative with respect to H a and H in the limit 
H a ,H b —> H, and identifying the singular contributions proportional to S(H a — H b ) on both sides of the equation 
lead to 

N^SlS{H-H a )=N f dSS 2 p(S,H) = -2G c d usp (k = 0;H). (44) 

The nonanalyticity in W[H, H; h] also implies that the 2-replica correlation function Gd,ij(H a , H a ; H b , H b ) contains 
a singular contribution proportional to S(H a — H b ). This induces singular contributions in the three disconnected 
direct correlation functions via the OZ equations and leads to formally diverging terms proportional to '(5(0)' when the 

replica fields are equal. The function Gd, however, has no obvious physical meaning and this problem is in principle 
harmless, although it may be a source of difficulties in an approximate treatment, as will be discussed below in section 
V. 

Finally, it is instructive to examine the behavior of the correlation functions in the limit H — ¥ ±oo, i.e. along the 
two branches of the hysteresis loop in the large-disorder regime. Assuming that the general behavior is the same as 
in the reference system (which seems quite reasonable, at least in the regime under consideration), we find 

G c = 0(e~^), C c = 0(H) 
G C = 0(1), C C = 0(1) 

G c = 0(H), C c = 0{e~\ & \) (45) 

and 

G d = 0{l), C d = 6(0)O(H 2 ) + O(H 2 ) 
G d = 0(H), C d = 0(H) 

&d = 5(0)O(H 2 ) + O(H 2 ), d d = 0(l), (46) 

where the notation indicates that both the regular and the singular contributions of Gd and Cd are of order H 2 . Note 
that we have considered the dependence on the source H. A similar dependence is found on rh when the latter goes 
to plus or minus infinity (recall that we work at the zeroth-order level of the expansion in free replica sums) . 

It is not surprising that G c vanishes as H — > ±oo. This is due to the already mentioned fact that there is only 
one metastable state along the loop at a given field so that the statistical fluctuations only come from the quenched 

disorder and are thus contained in the disconnected functions. (As a consequence, C c also vanishes.) The important 

feature is that G c and C c vanish exponentially fast with H. Indeed, this implies that the OZ equations for the 
physically relevant functions G c and Gd take in this limit the much simpler form: 

G c (k) = -J— 
C c (k) 

G d (k) = -S^-. (47) 

C c (k)2 

Note that Cdik) does not contain any diverging part when H — > ±oo (or m — > ±oo), as it must be. 
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IV. RANDOM PHASE APPROXIMATION (RPA) 



Equipped with the formalism of the previous section, we wish to introduce approximations based on assumptions 
on the structure of the direct correlation functions (proper vertices). Formally, these functions may be written as 



C = C (0) - s 



(48) 



where C^ ' is the 'bare' inverse propagator matrix (not taking into account the local potential V), whose components 
satisfy 



C^ )ab = 

M0)ab _ j, 
U ij — - J *iJ 

Cf ab = -A , 



Sab 



(49) 



and S is a 'self-energy' matrix (the minus sign is chosen so to match the usual definition of this quantity in field 
theory and many-body physics). In the above equations, Ay is 1 if i and j are nearest neighbors on the lattice, zero 
otherwise. 

From now on in this section, except if explicitly stated, we only consider uniform sources that are equal for all 
replicas. Our first approximation, akin to the usual random phase approximation (RPA), consists in assuming that 
the self-energies are identical to those of the noninteracting (J = 0) reference system. In Fourier space, this then 
leads to 



cr A (k) 


_ (jref 


C d (k) 


_ ^ref 
— °d 


cr A (k) 


= C c re/ -gJA(k), 


Q(k) 


_ firef 
- °d 


4 fliM (k) 


_ Are/ 
— c ' 


&(k) 


_ Are/ 



(50) 

where q = 2d is the connectivity of the hypercubic lattice and A(k) = (1/q) J^ e exp(ik.e) its characteristic function. 
The direct corre lation funct ions o f the refere nce system are obtained from the Green's functions computed in the 

appendix [Eqs. ( A18 |-( A20 ) and ( A27 |-( A37 | with K c = and K d = A]. The RPA effective action is immediately 
obtained by integrating the 'susceptibility sum-rule' d 2 (T / N) / dmdm — G~ 1 (k = 0) = C c (k = 0), which yields 

1 



^RPA 



N 



1 f 

(to, rh) — —T rel (m, m) — qJmrh. 



(51) 



From this expression, one obtains 



H(m, m) = H re f (m, rh) — qJm 
H (to, to) = H re f (m, to) — qjrh 



(52) 



It can be easily checked that T,Q PA (m, H)/N — ~T RPA (m 7 to) /N+rhH(m, to) is identical to the quenched complexity 
Eq (not to be confused with a self-energy) obtained in the mean- field model of Ref. [T^] (with qj replaced by J). 
(On the other hand, the auxiliary field H does not coincide with the parameter g introduced in this reference. In 
particular, H satisfies the Legendre relation 9(Eq(to, H)/N) /dm = —H, in contrast with g. In this respect, the RPA 
is a nicer way of obtaining the mean- field limit.) 

One expects the RPA to be valid when the coupling J is sufficiently weak, or, equivalently, when A and k are 
sufficiently large (all 'thermodynamic' quantities can be expressed in terms of the reduced variables V~A/J, k/J, and 
H/J). This is indeed what is observed in Fig. [2] where the predictions of the RPA for the magnetization along the 
ascending branch of the hysteresis loop are compared to numerical simulations (the descending branch is obtained 
using the symm etry H — > —H, to — > —m). The theoretical value is given by m RPA (H) = m re f (H re f) where m re f 

is given by Eq. (A21) with K c = and Kd — A, and from the above Eq. (52) the "displaced" field H re f is solution 
of the implicit equation H = H re ^ — q,Jm re f {H re f). One can see that the agreement is very good for J = 0.1 and 
deteriorates as J increases. In particular, the RPA overestimates the slope dm/dH and the magnetization curve 
already exhibits a reentrant behavior for J = 0.4 (typical of a mean-field theory below the critical point) whereas the 
actual system is still in the large-disorder regime. 
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FIG. 2: Magnetization along the ascending branch of the hysteresis loop: comparison of the predictions of the random-phase 
approximation (solid lines) to simulation data (circles) for various values of the coupling J (with k = 8 and A — 4). The 
simulation data are averaged over 1000 disorder realizations of linear size L — 30. 

The spin-spin and spin-random- field correlation functions along the loop for r = and r = e are shown in Figs. [3] 
and [I] They are denned by (see section III) 

G ss (r) ee sp7- m 2 = lim G d {r) (53) 

and 



Eqs. (47) and (50) then yield 



where 



and 



G sh (r) = s h r = A lim G c (r) . (54) 



Gf£ A (r) = GZ f P(r;z) (55) 
G PPA (r) = GL e/ [P(r; z) + zP'(r; z)\ (56) 

- = * RPA = (57) 



i r ... e ik - r 



p ^ = w?L/ v ^m (58) 

is the lattice Green function; in addition, P'(v;z) = dP(r; z) jdz (z and the correlation functions of the reference 
system are functions of m or H) . One can see that the RPA correctly predicts the values of the correlation functions 
along the hysteresis loop when J is small (we have checked that the agreement with the simulations remains good 
at the second and third nearest-neighbor distances) but considerably overestimates these values when J increases 
and the slope dm/dH becomes large (for instance around H 6.45 for J = 0.3). This can be traced back to an 
overestimation of the value of z (note that the RPA susceptibility diverges for z = 1). 

Finally, it is also interesting to examine the RPA correlation functions for H finite, which corresponds to metastable 
states inside the hysteresis loop in the H — m plane. For conciseness, we only consider the connected functions, as 
given by Eqs. [40j We then write 

c? PA (k)6? PA (k) - cf PA (k) 2 = c r c ef 6 r c ef - [c r c ef - <?JA(k))] 2 

= [C r c ef d r c ef - (C r c ef ) 2 ][l - *iA(k)][l - z 2 A(k)] (59) 
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H H 



FIG. 3: Spin-random field correlation function G s h(r) at r = (a) and r = e (b) along the ascending branch of the hysteresis 
loop for J = 0.1 (blue) and J = 0.3 (red). The predictions of the RPA (dashed lines) and the LSEA (solid lines) are compared 
to the simulation data (circles). 




FIG. 4: Same as Fig. 3 for the spin-spin correlation function G ss (r) 



where 



which yields 



- _ -RPA 
Zl,2 — 2=1,2 



cr f ± \/c c re/ 4 re/ 



Z\ 


zi - 


Zi 


zi - 


z 


z\ - 


Zi 


Z\ 




z\ - 


Zi 



Zi 


z\ - 


Zi 


zi - 


z 


z\ - 


~-2 


Zi 




z\ - 


~-2 



P(r;z 2 )] 

P(r;zi)] 

P(r;z 2 )] 



(60) 



(61) 



where z = z RPA as defined above. 
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It turns out that C r c e ^ is a positive function of H and H and G™' is negative, whereas the sign of C£ e ^ may change, 
as shown in Fig. [5] for H = 0. Therefore, depending on whether C r c e ^ is negative or positive, Z\ and z-i are real or 
complex conjugates, respectively. However, one can check that the correlation functions are always real, as it must 
be. On the other hand their behavior changes with the sign of C T c e ^ as illustrated by the leading asymptotic behavior 
as r — > 00. For instance, using the asymptotic expansion of P{n, n, n; z) |24j . one finds that 

Gf PA (n, n, n) ~ ^ \^^e~ x ^ - Il^± e ^ 
2ixn z\ — Z2 z\ — Z2 

(62) 

where Ai,2 = 31n[zi^ 2 /[l — Jl — zf i2 ] . Therefore, when C 7 c e f > and zi j2 are complex conjugates, the typical 

Ornstcin-Zcrnike fall-off e~^ n /n is modulated by oscillations at a wavevector Q = nQ^Ai^) where 9 denotes the 
imaginary part. For a given value of the magnetic field H, both the correlation length and the wavevector Q depend 
continuously on H (or, alternatively, on the magnetization m(H) of the metastable states), as shown in Fig. [6] This 
defines a 'disorder line' in the H — m plane where the qualitative behavior of the correlation function changes. In 
the case presented in Fig. [6] (for a small value of J for which the RPA is expected to be valid) , one finds that Q is 
non-zero for H = 0, that is for a typical metastable state at the field H . On the other hand, Q is always zero for 
H —> ±00, that is along the hysteresis loop. 





FIG. 6: Wave- vector Q as a function of the magnetization m(H) of the metastable states for J — 0.1 and H = 7. Q is non-zero 
for mtyp = m(H — 0), the magnetization of the typical states, but zero on the hysteresis loop (as H — > ±00). 
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V. BEYOND THE RPA 



A. Local self-energy approximation (LSEA) 

The RPA is exact in the limit of infinite dimension or infinite lattice coordination but, as we have just seen, is not 
a good approximation for d = 3 outside the weak-coupling regime. To go beyond this regime, we shall assume that 
the fluctuations renormalize the Green's functions without changing the functional form of their spatial dependence. 
This amounts to assuming that the corresponding self-energies are purely local. This assumption is the starting point 
of an approximate theory similar to the so-called Optimized Random Phase Approximation (ORPA) in liquid-state 
thcory[22 (see also 25J), the 'locator' approximation used in spin and Coulomb glasses [56- 28 , and to the dynamical 
mean-field theory (DMFT) for quantum problems [2 9], This type of approximation is known to be efficient when 
the physics is dominated by strong short-wavelength fluctuations rather than by long- wavelength fluctuations 29i. 
Accordingly, it cannot properly describe the vicinity of critical points but may capture for instance hysteresis and 
metastable effects that are already generated by the local potential in Eq. (J2j|. Note that the local self-energy 
approximation reduces to the RPA when d — > oo, but may be expected to provide reasonable results even in d = 3. 

From the point of view of thermodynamic consistency, it is useful to formulate the theory at the level of the 
effective action within the two-particle irreducible (2PI) formalism where T rep is considered as a functional of both 
the magnetizations {mf,rh^} and the Green's functions [30] (this is known as the 'entropy functional' in classical 
systems |31|). Up to an additive constant, T rep can be written as 

r rep [{m Q ,m a },GJ = S rep [{m a ,ih a }} + ^TrlnG^ 1 + ^TrCj°>£ + $,. ep [{m Q , m a }, g], (63) 



where G is defined in Eq. (36) and $ rep is the so-called Luttinger-Ward functional l?32l 133] . (<Iv ep is also called the 
2PI functional as it is in general the sum of all two-particle irreducible diagrams built with the fully dressed Green's 
functions. In the present case, the singular nature of the local potential makes the definition of the vertices appearing 



in the diagrams tricky and Eq. (63) must thus be considered as a nonperturbative definition of $ re p-) In the above 
expression, the trace involves a sum over both replica and spatial indices. The crucial point is that the explicit 
dependence of T rep on the pair interactions is contained in the classical action S rep [{m a , rh a }] and in the bare inverse 
propagator C' -'. On the other hand, for a given on-site potential, the functional $ rep is universal. By construction, 
the self-energies are obtained as 

d$ r 



= - 2 7^ ( 64 ) 



etc., and the stationnarity of the functional r rep against the variations of the Green's functions provides the 
Schwinger-Dyson equations, 

Q-l = £(0) _ J, (g 5 ) 



where one also has from Eq. (37) that G 1 = C. At the extremum, r rep then identifies with the physical (1PI) 



replicated effective action defineain Eq. (33) 



At the level of T rep , the local approximation consists in replacing Q rep by a sum of purely local contributions, 
Si 4>rep({ m i j in" }, G..). This readily implies that the self-energies are purely local. More specifically, restricting 
ourselves to the translationally invariant situation (uniform configurations), we find in Fourier space 

Q ab {k) = G (0)ab (k) - |T fc , (66) 
which expresses in a matrix form the equations for the components C ab , C ab , C ab in terms of their counterparts 



Q(o)ab an( j j n k e self-energy. According to Eq. (49), the only nonzero components of C^ ab are C^ ab (k) = 

-qJ\{k)S ab and C^ ab {k) = - A. The self-energies are then functions of {m a , m a }, J, k and A to be determined. 

As discussed for instance in Ref . [35] , the simplest strategy for computing the above functions is to define a single-site 
effective action (the so-called 'impurity' model in strongly correlated Fermi systems) that involves the original on-site 
interaction and arbitrary quadratic terms. This model is in general exactly solvable and self-consistency equations are 
then obtained by imposing that the Green's functions of the single-site action coincide with the site-diagonal Green's 
functions of the original lattice model. 

We thus introduce the single-site action 

S a ,re P [{s a ,s a }] = -\j2 [K ab s a s b + K ab {s a s b + s a s b ) + A ab s a s b ^ +^s'V(s a ) (67) 

a, fa a 
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and add external sources Hq , irfg that fix the same 'magnetizations' rh a and m a , respectively, as in the fully interacting 
system. The so-called 'Weiss fields' K ab ,K ab ,K ab must be chosen so that the Green's functions Gff,6tf,Gf of the 



effective single-site model coincide with the on-site Green's functions Gf b , Gf b , Gf b of the original model, i.e. 



G ab (r = 0) = G 



ah 



(68) 



with identical self-energies Ef b = E ob %, T,f b = £ ob <%, £? b = E ab % [34] . By definition of the self-energies, the direct 
correlation functions of the effective model are given by 



(69) 



as = i n the single-site model, with the matrix K ab collecting the Weiss fields K ab , K ab , K ab . Subtracting 

from Eq. ( 66 ) then yields 

C Qb (k) = C (0)ab (k) + K ab + C a Q b 

and the consistency requirement between G° b and G ab provides the following equation, 

-l 



(70) 



G 



= o 
i-i 



rfk 



C (0) (k) + K + C 



(71) 



where frfk is a short-hand notation for f d 3 k/(2ir) 3 . 

Eq. ( 71 1 can be considered as a self-consistent equation for K ab since the direct correlation functions of the single- 
site effective model are obtained from the action in Eq. (67) as functions of the replica 'magnetizations' and of the 
Weiss fields: C ab = C ab [{m e ,rh e }; {K e ^}]. Again, we can use the expansion in number of free replica sums to derive 
explicit expressions. We decompose the matrices according to K ab = K°<5 & + K^ b , etc., and expand each component 
as m Eqs. |31p2| (fi xing the magnetizations instead of the sources). At zeroth-order, one finds from Eq. (|7l|) 



G n (m Q ,m a ) = C n (m a ,m a )- 1 = dk 



C (0) (k) + K (m a ,m a ) + C (m a ,m a ) 



(72) 



and 



-g 0(j (m", m a ; to , m b ) = Q Q c (m Q , m a )~ Q Q d {m a , m a ; m b , m b )Q Qc (m b , m b ) 



b\-l 



rfk 



C (0) (k) + K (m a ,m a ) + C n (m a ,m a ) 

1 



-(0) 



Kf a - a b ~ b\ , /~i / a ~ a b ; b\ 
d (m , m ;m , m ) + C d {m , to ; to , to ) (73) 



G^(k) + K C (™ > b ) +Q 0c {m b ,m b 



where C c (m a , rh a ) and C rf (w a , m a ; to , to ) are calculated from the single-site effective action in Eq. (67) with the 
Weiss fields now functions of the magnetizations through the above implicit equations; the latter are considered at 
zeroth-order in the expansion in number of free replica sums, i.e. K° = K Jjn a , rh a ) and K"^ = K d (m a , to q ; TO b , TO b ). 

As the Luttinger-Ward functional is universal, the expression of T rep is obtained by subtracting from Eq. ( 63 1 the 
formal expression of r rePi o for the single-site effective model. Restricting again the calculation to the translationally 
invariant situation (uniform configurations), we obtain 

j^T rep = T repfi - qj^m a rh a + \Y1 [K ab m a m b + K ab (m a m b + m a m b ) + [k ab - A)m a m b 

a a,b 

- ^ y rfk [lHn (G(k)Go 1 ) - Tr G (0) (k)G(k) + Tr G ( 0) G ] ■ (74) 

This effective action, when evaluated at the extremum corresponding to the solution of the Schwinger-Dyson equations, 
is a function of the replica magnetizations that it can also be expanded in number of free replica sums, 



r rep ({m a , m a }) = ri(TO Q , to q ) - r 2(™ a ' ™ Q ; 



TO b ,TO b ) + 



(75) 
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the pth term then providing information on the pth cumulant (see above) . 

Since the local approximation is performed at the level of the Luttinger-Ward functional, the theory possesses some 
form of thermodynamic consistency[3S]. Indeed, one obtains the same value of T rep by using the above equation or by 
integrating once a quantity that depends on the Green's functions (e.g. the 'internal energy') with respect to the model 
parameters (e.g. the coupling constants). For instance, using the fact that the functional is extremal with respect 
to the variations of the Green's functions and considering for simplicity the case where all replica magnetizations are 
equal, one readily derives the relations 



dT 1 /N _ 

dYx/N _ 1 
<9A ~~ ~2 



G(r = e) + mi 

d c (r = 0), (7(»i 



which are exact and are thus preserved by the local approximation. The extremal property also allows us to simply 
relate the fields H(m,rh) and H(m,m) to the fields Hariri, rh) and HQ(m,rh) in the single-site effective system (they 
all correspond to zeroth-order terms in the expansion in free replica sums) by only considering the explicit dependence 
of F and Tq on m and rh. This yields 



H(m, rh) 



1 _ 1 dS(m,m) 
N drh N dm 



tt , ^_^ r o,ii dS (m,m) 

H {m,m) = — '- \ = — , (77) 

am 1 Ik am 

which leads to 

H(m, rh) = Ho(m, rh) + [K c (m, rh) — qj]m + K c (m, m)m . (78) 

Similarly, one finds that 

H(m, rh) = Ho(m, rh) + [K c (m, rh) — qj]rh + K c (m, m)m . (79) 

One may however notice that the local app roximatio n does not really yield a fully consistent theory. For instance, 
the value of C c (k = 0) obtained from Eqs. |66|69|71 1 does not coincide with that obtained from the 'susceptibility 



sum-rule' C c (k = 0) = d 2 (T i / N) / dmdrh — dH(m,m)/dm. This inconsistency is a well-known flaw of this type of 
approximation and it can be cured by renormalizing the value of the direct correlation functions at nearest-neighbor 
distance as is done for instance in the so-called Self-Consistent Ornstein-Zernike Approximation (SCOZA)|251 136] . 
This route, however, looks prohibitively difficult in the present case and will not be pursued. Improving the theory by 
introducing a k-dependence in the self-energies as is done in the cluster dynamical mean-field thcory[29 can make the 
theory exact up to order 1 /d, which is an interesting property, but it does not solve the above inconsistency problem. 
At this stage, we must point out a serious difficulty occuring in the approximate framework developed above. We 

have stressed that, just like the reference system, the exact system is most likely such that the Green's function Gd 
is singular with a diverging term proportional to 5(0). The 2PI functional and the other Green's functions being 

most likely finite, the contribution of Gd must altoget her vanish in the expressions of these finite quantities. One 



can therefore plainly drop all dependence on Gd in Eq. (63), which amounts to a perfect cancellation between terms 

in TrlnG -1 and in $ rep (there are no Gd contributions in TrC^G). This of course has consequences on the self- 
energies. To avoid inconsistencies, this property must be satisfied, or at least enforced, in any sensible approximation. 

As shown in the appendix, it is easily realized that Go,d in the single-site model may develop a singular £(0) term 

only if the Weiss field K c is identically zero. A sensible approximation scheme must therefore be compatible with (i) 

setting K c to zero, (ii) discarding the self- consistent equation on Gd, and (iii) dropping all contributions involving 

Gd in the two-particle irreducible functional of the original model and of the effective single-site model. It turns out 
that these requirements are not met by the local self-energy approximation inside the hysteresis loop (when H and 
m are finite and the quenched complexity is strictly positive). Along the hysteresis loop (H, rh —¥ ±oo), the situation 
improves and the local self-energy approximation is well behaved, at least at the level of single-replica quantities. 
Awaiting for a proper resolution in the general case (see the discussion in conclusion), we now discuss the behavior 
on the hysteresis loop. 
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B. Behavior along the hysteresis loop 



We consider the solution of the LSEA discussed above along the hysteresis loop, when H — > ±00 or alternatively 
when m — > ±00. Our main assumptions, motivated by the expected exact behavior along the hysteresis loop, are 

that the Weiss field K c and the corresponding direct correlation function Cq jC in the effective model are identically 
zero (see above) and that the complexity is zero as well. The former of these two assumptions implies from Eq. (40) 

that G c (k) and thus G c (k) in the original model are also zero, as anticipated. We are then only interested in the 
correlation functions that are physically observable, namely G c (k) and G^k). As a result, we only need to solve 
the two corresponding self-consistency equations, G c (r = 0;m a ) = Go jC (m a ) and G d (r = 0;m a :m b ) = Go t d(rn a ;m b ), 
leading to 



G , c (m a ) = / dk 



1 



G , c (m Q ) + K c (m a ) - q JA(k) 



(80) 



G 04 {m a ;m b ) = - / dk 



G M (m a ; m b ) + K d (m a ;m b ) - A 

'[CoA™ a ) + K c (m a ) - qJ\(k)}[CoA™ b ) + K c {m b ) - qJX(k)} 



(81) 



where Go, c (w), Go, c (m) = G l(m), Go, d (m a , m b ), and Go,d(r7j a ; m b ) = — Go, c (™°)Go,d("i a , m b )Co tC ( mb ) are obtained 
from the effective single-site action 



S ,re P [{s a , s a }} = [-Ks a + V'(s a )] s a - ^ kfs* 



(82) 



with the magnetizations m a fixed by the external sources Ho(m a ) [rh a — > ±00 follows from Ho(m a ) — > ±00) and the 

Weiss fields implicitly determined as functions of the magnetizations, = K c (m a ), Kf 1 = K,i(m a ;m b ). The actual 
computation of the correlation functions associated to this effective single-site action is performed in the appendix. 
After introducing 



z(m) 



qj 



G , c (m) + K c (m) 



P{z) = P(r = 0; z), and P'(z) = dP{z)/dz (as in section IV), Eqs. ( |80| and (81 1 can be rewritten as 

1 _ 1 

G , c (m) G , c (m) + K c (m) 

where we have dropped the superscript a on the magnetization, and 



-P(z(m)), 



C 0id (m a ;m b ) 



Go, d (r7i a ; m b ) + K d {m a ;m b ) - A 



G ,c(m")Go, c (m fc ) [G , c (m«) + X c (m-)][G , c (™ b ) + K c (m b )} 
For equal magnetizations m a = m b — m, Eq. (85) simplifies to 



z{m a )P{z{m a )) - z{m b )P{z{m h )) 
z(m a ) — z{m b ) 



(83) 



(84) 



(85) 



C 0)rf (m;m) G ,d(m; m) + K d (m; m) - A 



G , c (m) ; 



[Go, c (m)+X c (m)] 2 



[P(z(m)) + z(m)P'(z(m))] . 



(86) 



Once the coupled equations (84) and (86) are solved for K c (m) and K d (m;m), we calculate the field H along the 
loop from Eq. (78 1 with K c = 0, 



H(m) = H (m) + [K c (m) - qJ]m, 



(87) 



and we obtain the physical correlation (Green's) functions for equal magnetizations from Eqs. (47) as 

P(r;z(m)) 



G sh (r;m) = A G , c (ra) ; 



P(z(m)) 
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FIG. 7: Magnetization along the ascending branch of the hysteresis loop: comparison of the predictions of the local self-energy 
approximation (solid lines) to the simulation data (circles) for various values of the coupling J. The dashed lines for J = 0.4 
and J = 0.5 indicate the macroscopic jump predicted by the theory. 



G ss (r;ra) = Go,d{m;m) 



P(r; z(m)) + z(m)P'(r; z(m)) 
P(z(m)) + z(m)P'(z(m)) 



(89) 



with explicit expressions for Go, c ( TO ) & n d Grj,d(TO;rn) given in the appendix . 

The predictions of Eq. (87 1 for the ascending branch of the hysteresis loop are shown in Fig. [7] Comparing with 
Fig. [2j we see that the improvement over the RPA is quite significant: the agreement with the simulations is now 
satisfactory up to J = 0.4 where a small reentrant behavior is observed in the upper part of the curve. This behavior, 
which erroneously indicates that the system has entered the small-disorder regime, is related to the fact that the 
effective action is obtained via the 'energy' route which is not consistent with the 'susceptibility' one, as mentioned 
above. This lack to thermodynamic consistency as the coupling increases is illustrated in Fig. [HJ Of course, the actual 
magnetization cannot decrease as H increases and it must jump at a 'spinodal' field where the slope dm/dh diverges 
for the first time (as can be seen in Fig. [7]for J = 0.5, the theoretical curve may display several spinodal fields). Note 
that due to this inconsistency the condition dm/dH — > oo does not imply G c (k = 0) = and thus not z = 1 (from 
Eq. S one has G c (k = 0) = (l/A)G sh (k = 0) = G 0>C /[P(»(1 - z)}). 




FIG. 8: Test of thermodynamic self-consistency in the LSEA: comparison of the susceptibility G c (k = 0) (solid lines) with the 
slope dm/dH of the magnetization curve (dashed lines) along the ascending branch of the hysteresis curve. For J = 0.4 the 
slope diverges at the spinodal field H sp m 6.08 whereas G c (k = 0) remains finite. 
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The correlation functions G s / l (r) and G ss (r) at r = and r = e are compared to the simulation data and to the 
RPA in Figs. [3] and [3] The agreement with the simulations is now excellent for J = 0.3. In particular, the LSEA 
correctly predicts that the maximum of G s f l (r=0) decreases as J increases, contrary to the RPA. The values at the 
second and third nearest-neighbor distances[37] are also shown in Figs. [9] and [To] The theoretical predictions for 
J = 0.3 are still fairly good although the values of the functions are slightly overestimated. More generally, the results 
displayed in these figures show that the actual dependence of G s ? l (r) and G ss (r) with distance is well described by 
RPA-like expressions in the large-disorder regime and not too close to criticality. The two functions involve a single 
correlation length £ [e.g. the second-moment correlation length defined by G s ^(k) ~ G s h(0)(l + £ 2 fc 2 ), k — > 0, and 
thus related to z by qt; 2 = z/{\ — z)] and the main effect of disorder fluctuations is to renormalize £. 




FIG. 9: Spin-random field correlation function at the second (a) and third (b) nearest-neighbors along the ascending branch 
of the hysteresis curve for J = 0.1 (blue) and J = 0.3 (red). The predictions of the LSEA (solid lines) are compared to the 
simulation data (circles). 




FIG. 10: Same as Fig. [9] for the spin-spin correlation function G ss (r). 

From the above correlation functions, we can compute the average energy per spin along the loop, which is the sum 
of four contributions, 



U/N = H/N = mi + u 2 + u 3 + u A 



(90) 
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with 



Ul 

«4 = 



J_ 

N 



E 



qj 



<hJ> 
—~s~iH = —mH 

-7JTi = -G sh (r = 0) 



V(s t ) = ±-(H + JJ2 S J + h if = ^i H2 + A + ZqJmH + 2qJG sh (r = 1) 



j/i 



2k 



+ qJ 2 (G ss (r = 0) + to 2 ) + q(q - 2) J 2 (G ss (r = y/2) + to 2 ) + qJ 2 (G ss (r = 2) + to 2 )] 



(91) 



where, again, the dependence of H and of the Green's functions on the magnetization to is left implicit. As can be 
seen in Fig. [TTJ the average energy per spin is also very well reproduced for J = 0.3 and the discrepancies for J = 0.4 
are limited to the small range of H where the reentrant behavior in the magnetization occurs. 




10 12 



FIG. 11: Energy per spin along the ascending branch of the hysteresis loop. The predictions of the LSEA (solid lines) are 
compared to the simulation data (circles). 



C. Two-replica correlation function and avalanches along the hysteresis loop 



Additional information on the system along the hysteresis loop is encoded in the 2-replica spin-spin correlation 
function for distinct magnetizations, G ss (r;m a ,m b ) = lim - •'*•• •••"•» — <""'■• < '-"< ■ rr'- 

From the equations above, one derives 



G ss (r;m a ;m b ) = G , d (m a ;m b ) 



H^±o G d(v;m a ,m b )) = s (H-)s r {H b ) - m{H a )m{H b ). 
z{m a )P{v-z{m a )) - z{m b )P{v-z{m b ))~ 



z(m a )P{z(m a )) - z(m b )P(z(m b )) 



(92) 



which gives back Eq. (89) when the magnetizations are equal. As discussed in section III, this function allows 
one to compute the unnormalized second moment of the avalanche-size distribution via Eq. ( 44 1 , provided one 
extracts the nonanalytic cusp-like dependence of G ss (k = 0, m a ; m b ) on the difference m a — m b , and consequently on 
H a — H b . The amplitude of the linear-cusp contribution is then obtained from dG ss (k = 0;m a ;m b ))/d(m a — to 6 ) 
when (to" — to 6 ) — > + . It is clear from Eq. (92) that the cusp only comes from Go } d(iri a ; to 6 ), i.e. from the two- 



replica correlation function in the single-site effective model defined by the action in Eq. ( 82 ) . This function depends 
on the magnetizations both in an explicit way and through the Weiss fields, which we write as Go j£ i(TO a ; to 6 ) = 

Go : d(TO a ;TO 6 ;i^ c (m Q ),i^ c (TO 6 ),iY d (m a ;TO a ),i^ d (TO 6 ;TO 6 ),/^(TO a ;TO 6 )). The contribution of the Weiss fields to the 
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cusp can only come from K d (m a ; m b ) [with m a ^ m b ], implying that 



dG „ 



d(m a 



dGo.d 



o+ 



d(m a — m b ) 



k a ,k dl o+ dK d (m a ;m b ) 



dK d (m a ;m b ) 
o+ d(m a 



m b ) 



(93) 



o+ 



where the subscript + indicates that the derivatives are evaluated for (m a — m b ) — > + and the derivative with respect 

to K d (m a ; m b ) in the right-hand side is taken with m a , m , K c (m a ), K c (m b ), and K d {m a ; m a ), K d (m b ; m b ) fixed. On 
the other hand, the self-consistency equation (81) imposes that the cusp in Go,d(m a ; m b ) is directly proportional to 

the cusp in the Weiss field K d (m a ;m b ). Indeed, replacing Co, d (m a ;m b ) by —Co, c (m a )Go, d (m a ;m b )Co, c (m b ) in the 



right-hand side of Eq. (81 1, one finds after simple manipulations 



dG ,d 



o+ 



G 0td (m,m) 
K d (m; m) - A. 



dK d {m a ;m b ) 



d(m a — m b ) 



(94) 



o+ 



d(m a — n 

where m is the common limit of m a and m b . 

Before we proceed any further, we must however deal with a difficulty that arises in Eq. (93) from the behavior 

of the partial derivative dGo, d /dK d (m a ; m h ) when (m a — m b ) — » + . At zeroth-order in the expansion in number 
of free replica sums, it is equivalent to use the replica magnetizations ({m a }, {m a }) or the corresponding sources 
({Hq = Ho(m a )}, {Hq = Ho(m a )}), so that we may as well consider the partial derivative of Gq^H^Hq) with 

respect to K d (HQ-, Hq). This quantity is most easily expressed by using the property that Go, d (H§ , H§ ; Hq, Hq) is 
the second derivative of the second cumulant W 2 (Hq , Hq ; Hq, Hq) with respect to Hq and Hq. On the other hand, 
from the definition of the action of the effective model, the derivative of W2 with respect to Kf — Kd{HQ,HQ) is 

equal to & 04 (H%,H%] H b , H b ) + m(H%, H^)fh(H b , H b ) (using the symmetry kf = K ba ). After permuting the order 
of the derivatives, one thus has 



dGg : d 
dk d (Hg;H b ) 



d 2 

8H^dH b 



[G 0td (HS,H$; H b , H b ) + m(F °, H-)m(H b , H b )} 



(95) 



is 



a quantity that is still defined i n the limit Hg Hft ±00 since G , d ( H o > H o > H o> H o) + m ( H o> ^0)^(^0) H o 
proportional to H§Hq [see Eqs. (A33) and (A34)]. The problem is that the Green function Go,d contains a singular 
term proportional to S(Hq — Hq), as already mentioned, so that dGo y d/dK d (H^; Hq) diverges when Hq — H b — > 



or, equivalently, when m a — m — > 0. This causes the right-hand side of Eq. (|93j) to diverge, whereas Eq. (94) does 
not predict any divergence. This spurious behavior is intrinsic to the local self-energy approxi mat ion but it can be 

circumvented by simply discarding the singular contribution of Go,d & n d only keeping in Eq. ( 95b the contribution 
of the regular term. This is admittedly a drastic regularization procedure which however may be rationalized by 

recalling that the singular term, which is an exact feature of the Green function G d , should not play any role in an 
exact treatment [4*0]. Then, combining Eqs. (93) and (94) and considering the fields instead of the magnetizations lead 
to 



dG Q . d (Hg ; Hq) 



d(H- - H b ) 



0+ 



k d {H ;H )-A 
Go,d(Ho, H ) 



dGo,d 



dK d (H*;H b ) 



dGo. d 



0+ 



H b a ) 



(96) 



K a ,K d ;0+ 



where Hq = H (m) is the common limit of Hq and Hq. It is important to note that the two derivatives of Gq^ that 
appear in the right-hand side can now be computed from the effective single-site action in Eq. ( 82 ) by assuming that 



the Weiss fields are replica-symmetric and independent of the sources Hq and Hq. This calculation is performed in 
the appendix. 

Once the amplitude of the cusp in \m a — m b \ is known, the amplitude of the cusp in \H a — H b \, G™ sp (k — 0,H), 
which is needed for computing the unnormalized second moment of the avalanche distribution via Eq. ( 44 ) , is obtained 
through 



GT s sp (k = 0,H) = 



dG ss (k = 0;m a ;m b ) 



m b ) 



0+ 



dm 
dH 



(97) 
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with H and m related via Eq. (87). One finally obtains from the Fourier transform of Eq. (89 1 

dG Q , d (m a ;m b 



G 



cus P(k = 0,H) = 



d(m a — m b ) 



d(H$ - H b ) 



(dm/dH) 

0+ [1 - z(m)) 2 [P(z(m)) + z(m)P'(z(m))] 

(dm/dH)/(dm/dH ) 

0+ [1 - z(m)] 2 [P(z(m)) + z(m)P'(z(m))} ' 



(98) 



with dGo,d(HQ ; Hq)/8(Hq — Hq)\ + given by Eq. (96 1. This equation further simplifies if one replaces (dm/dH) by 
G c (k = 0) = G , c /[(1 — z )P( z )]j thereby neglecting the small thermodynamic inconsistency of the present approach. 
(On the other hand, (dm/dH ) is exactly equal to Gq, c in the single-site effective model.) 





2 4 6 8 10 12 2 4 6 8 10 12 

H H 



FIG. 12: The unnormalized second moment of the avalanche-size distribution along the ascending branch of the hysteresis 
curve for J = 0.1 (a) and J = 0.3 (b) . The predictions of the LSEA (solid lines) are compared to the simulation data (circles) 
resulting from an average over 500 disorder realizations of linear size L = 10 for J = 0.1 and L = 20 for L = 0.3 |16j . 

The theoretical prediction for the unnormalized second moment of the avalanche-size distribution is compared to 
the simulation data in Fig. |12| One observes a reasonable agreement, even for J = 0.3, which gives some a posteriori 
justification to the regularization procedure that has been used in the above calculation. 



VI. CONCLUSION AND OUTLOOK 



In this paper we have proposed a formalism to study the out-of-equilibrium hysteresis behavior of random field 
systems when submitted to an adiabatically varying external field at zero temperature. The key ingredients consist 
in relating the out-of-equilibrium behavior to the statistics of the metastable states, introducing auxiliary variables to 
handle the latter, and building an approximation scheme based on the structure of the pair correlation functions. We 
have applied this program to a soft-spin version of the random field Ising model and focused on describing the system 
along the hysteresis loop, which is the envelope of the metastable states in the magnetization-applied field plane. 
The physics of the problem involves 'avalanches' between metastable states that in turn generates nonanalyticities 
in the dependence of several correlation functions on their arguments. We have used an approximation borrowed 
from condensed-matter theory and most conveniently formulated in a 2PI framework, whose lowest order amounts 
to neglecting the spatial dependence of the self-energies. We have derived in this way predictions for the physical, 
spin-spin and spin-random-field, pair correlation functions, as well as for the second moment of the avalanches, 
which can be compared to computer simulation data: we find a good agreement between predictions and simulation 
data above the (out-of-equilibrium) critical point, which represents a significant improvement over the mean-field 
(random-phase approximation) results also calculated here. Away from criticality, the correlation functions along 
the hysteresis loop keep essentially the same spatial structure as in the random-phase approximation with however a 
strong renormalization of the correlation length due to disorder-induced fluctuations. 



22 



In spite of the encouraging accuracy of the results, we have encountered difficulties in extending the basic local 
self-energy approximation either to the situation 'inside' the hysteresis loop, for which the complexity associated with 
the number of metastable states is nonzero, or to the many-replica correlation functions, even along the hysteresis 
loop. The presence of avalanches at T — indeed generates a strong singularity in the form of a Dirac delta function 
in the (unphysical) correlation function of the auxiliary variable. The contribution of the latter must therefore exactly 
vanish in physical quantities and in the effective action. Such exact cancellations are however hard to implement in 
an approximate theory. To cure this problem, one must somewhat relax the self-consistency of the local self-energy 
approximation to allow for one more constraint enforcing the vanishing of the diverging terms without overconstraining 
the theory. In addition, to improve the accuracy of the predictions and, for instance, to provide a good description 
of the much studied hard-spin Ising model in a random field, one will have to go beyond the local approximation 
for the self-energies, as done in the cluster dynamical mean-field theory Work in this direction is in progress. 
The vicinity of the out-of-equilibrium critical point on the other hand cannot be treated in such cluster extensions of 
the local approximation and requires a different treatment along the lines of the recently developed nonperturbative 
functional renormalization group [20 . 



Appendix A: Single-site effective model with 'replica-symmetric' Weiss fields 

In this appendix we compute the correlation (Green's) functions of the single-site effective model whose partition 
function in replica space reads 

Z rep ({H a ,H a }) = f ]Jds a ds a e-S-pKsV^+EJ^+^s ] (Al) 

where the action is given Eq. ( 67 ) (hereafter, for ease of notation we drop the subscript on all quantities.) We consider 
a replica-symmetry ansatz for the Weiss fields, K^ h = K 8nh+K with K_ and K. taken as fixed. Although this ansatz 
neglects the fact that the Weiss fields depend on the magnetizations {m a } and {rh a } through the self-consistency 
equations ( 72 1-( 73 ) , the results derived in this appendix are nonetheless sufficient to compute all the quantities studied 
in the present work, including the amplitude of the linear cusp in the two-replica spin-spin correlation function. The 
action then becomes 

S rep ({s a , s a }) =-\j2 [ K ^ S<1 ) 2 + 2K c s a S a + K c (s a ) 2 } - \ [K d u 2 + 2K d uv + K d v 2 ] + s a V(s a ) (A2) 

a a 

where u = ^2 a s a and v — ^2 a s a . The quadratic dependence on u and v can be eliminated by using a Hubbard- 
Stratonovich transformation with two auxililary fields hi and h 2 that play the role of correlated random fields. This 
yields 



Z rep ({H a ,H a }) = 1 -= f dh 1 dh 2 e-^ kdhl ~ 2kdhlh2+Kd ' 1 ^ TT / ds a ds a 

2tt v R d J A A J 



where 



f a (s a , s a ; h u h 2 ) = X - [K c {s a f + 2K c s a s a + R c {s a ) 2 ] + s a (H a + h x ) + s a [H a + h 2 - V'{s a )] (A4) 



and R d = K d K d — K\ must be a positive quantity. 

When all sources act identically in each replica (i.e. H a = H, H a = H), we then have 

Wi(H,H) = lim - ]nZ rep (H,H) = — 1 -= [ dh x dh 2 e~^ Kdhl ^ Kdhlh ^ Kdh ^^{B,B\h x M) (A5) 
«->-o n 2ir\fR d J 

with 

W(H, H; hi, h 2 ) = \njds ds e f{s ' S ' MM) . (A6) 

To compute this quantity we first integrate e-^ s ' s ) over s along the imaginary axis (taking into account the factor 
l/(2z7r) that was adsorbed in ds). This gives 

ds e f(s,s-MM) = _1_ exp ( W!> ) (A7 ) 
2nK r 2 ^c 



23 



with 

g(s; hi,h 2 ) = K c K c s 2 + 2K C (H + hi)s -[H + h 2 + K c s - V '(s)} 2 

= [K C K C - (k - K c ) 2 ]s 2 + 2[(H + hi)K c + (k- K C ){H + h 2 + k sign(s))]s - (H + h 2 + k sign(s)) 2 . 

(A8) 

The second integration over s from — oo to and from to +oo finally yields 

W(H, H; h u h 2 ) = -\ ln(jy + In C ~ "f^ e x ~ + 1 ± °f Y+) e*+) (A9) 

with 

X± = ^(H + h 1 ) 2 + - i -(H + h 2 ±k)lK c {H + h 2 ±k)+2(k-k c )(H + h 1 )} (A10) 

and 

y± = , ] [k c (k + in) + (k- k c )(h + h 2 ± k)} . (ah) 

2K C R C 

Like Rd, the quantity R c = (k — k c ) 2 — K c k c must be positive. The Green's functions can be calculated by derivation 
of Wi(H, H) with respect to the Weiss fields, 

— = -{G c + G d + m>), ^ r -f c 

— — = G c + Grf + mm , — — = G c 
dK c dK d 

= o(G c + G d + m 2 ) , — ^- = -G c , (A12) 
2 2 



and the corresponding direct correlation functions are then obtained by inverting the Ornstein-Zernike equations ( 40 ) 



and (411 



It is rather obvious that this set of equations leads in general to an analytic behavior of the correlation functions as 

a function of the source H . In particular, the function Gd does not have a singular term proportional to (5(0) although 
this is the expected behavior in the original lattice model, as p o inted out in the main text. One can easily see that 



the condition for a singular behavior to emerge from Eqs. (A9)-(All ) is that K c = 0: the quantity Y±(H, H; hi, h 2 ) 
then goes to ±oo depending on the sign of H + h 2 ± k and this induces a Heaviside step function in Eq. ( A9 1 and a 
Dirac delta when differentiating with respect to H. 

We now focus on the behavior along the hysteresis loop where it is sufficient to only keep the two Weiss fields k c 

and Kd (however, the limit H — > — oo which corresponds to the ascending branch will only be taken at the end of the 
calculation). The starting point is the simpler partition function 

Z rep {{H a },H a }) = f Y[ds a ds a e^£".* lV e T, a {s a H a +s a [Kc3 a +H->-v'(s°-)]} ( A13 ) 

which also encompasses the case of the reference system (J = 0) where k c = and k d = A. A single auxiliary 
(random) field h is now required to eliminate the quadratic dependence onv = ^ a s a . This leads to 



Wi(H,H) = J dhp(h)W(H,H;h) (A14) 
where p(h) is a Gaussian distribution with zero mean and variance Kd (which plays the role of a 'renormalized' 
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disorder) and 

W(H, H- h)=\nj ds ds e sH+'s[k cS +H+h-v'(s)] = i n J ds e sH § ^ k _ x c ) s -H-h-k sgns] 

= - ln(fc - K c ) + H H + H 7 k + [-^- + In (2 cosh -^L-)]6(H + h + k) 
k - K c k - K c k — K c 

kH kH 

+ [ ^- In (2 cosh ^)]e(H + h-k) . (A15) 

k — K c k - K c 

As a consequence, we find 

- H — k 9k H 
W l {H, H) = - ln(fc - K c ) + H £ + - ^V(H + k) 

k - K c k — K c 

k kff kH 

+ , -[In (2 cosh ^) - [P(H + k)- V(H k)} (A16) 

k - K c k — K c k — K c 

where V{x) = p(y)dy — J^° x p(y)dy = |[1 + erf— t==]. From this, we readily obtain the magnetizations 

V 2K d 

f)Wi H — k 9k k kH 

m(H, H) = °^ = — y + —^-V(H + k) + — -[tanh —- 1] [P(H + k) V{H k)\ 
oH k — K c k — K c fc — K c k — K c 

m (H,H) = ^ = — * + ^Lp{H + k) + [In (2 cosh -™) -^L] [p(H + k) p(H k)\ , (A17) 

On k — K c k — K c k — K c k — K c 

and the connected (Green's) correlation functions 

G C (H, H) = ^= - ' [1 tanh 2 -^][P{H + k) - V(H k)] (A18) 
an (k — K c y k - K c 



f>rh 1 9k k kH 

G C (H, H) = - s = + — - ¥ p{H + k) + [tanh — - y 1] ]p(H + k) p(H - k)} (A19) 

an k — K c k — K c k — K c k — K c 

s, - dm 2kH H + k kH kH H + k H-k 

G C (H,H) = — = p(H + k)- [ln(2cosh - tt\[— — p(H + k) — p(H - k)] . 

oH (k-K c ) Kd k-K c k-K c Kd Kd 

(A20) 

In the limit H — > — oo, these expressions simplify to 

m{H) = — l —^[H -k + 2kV{H - k)] 
k- K c 

m(H,H)~ — r ^ r [l + 2kp(H-k)} , (A21) 

k - K c 

and 



G C (H) = 

G c {H) = —^[l + 2kp{H-k)\ 
k-K c 

s - 2kH 

G c {H,H)~-„ (H-k)p(H-k). (A22) 

K d (k - K c ) 

At zeroth-order of the expansion in number of free replica sums, the disconnected functions are obtained 
from the derivatives of the second cumulant W 2 (H a , H a ; H b , H b ) = J dh p(h)W(H a , H a ; h)W{H b , H b ; h) - 
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W 1 [H a ,H a ] W 1 [H b ,H b ] with respect to the sources. We first consider G d {H a , H a ; H b , H a ) = 
d 2 W 2 (H a ,H a ;H b ,H b )/dH a dH b and 

G d (H«,H«;H b ,H>>) + rn(H",H«)rn(H b ,H b ) = [ dhp(h) dW{H ^ h) ' k) ■ (A23) 

J aH a oH" 

The presence of Heaviside step functions in W(H,H;h) makes the result dependent on the sign of H a — H b . For 
simplicity, we set at once H a = H b = H. We then find 

G d (H a , H a ; H b , H b ) + m(H a , H a )m(H b , H b ) 

= K [f sym (H a , H; H b , H) + g{H a , H)6(H b - H a ) + g(H b ', H)6(H a - H b )] 

(k — K c ) 2 

where f sym is a symmetric function of H a and H b , 

f s v m (H a , H; H b , H) = K d + (H a - k)(H b - k) 

kH 

+ k[l + tanh H \(H b - k)V{H a + k) + (H a - k)V{H b + k) + K d [p(H a + k) + p(H b + k)}] 

k - K c 

kH 

+ k[l - tanh — ] [(H b - k)V(H a - k) + (H a - k)V(H b - k) + K d [p(H a - k) + p(H b - k)}] 

k — K c 

kff 

+ k 2 [l- tanh 2 ^][V(H a -k)+ V{H b - k)] , (A25) 

k K r 



and 



kH kH 

g(H, H) = k 2 [l + tanh ^-] 2 V{H + k) + k 2 [l - tanh ^] 2 T(H - k) . (A26) 

fc — K c k — K c 



Therefore, when H a ,H b — >• H, we obtain 



1 



G d (H a , H; H b , H) + m(H a , H)m(H b , H) = - \ H; H, H) + g(H, H) 

(k - K c y 



2 1 a 0| dH 

For H a = H b = H -> -co, G d (H; H) is then given by 



- 1 -\H a -H b \ d9 ^ H l+0((H^-H b ) 2 )] . (A27) 



G d {H; H) = Rd[1 + 4kp(H k)] + 4k2r(H k)[1 V{H k)] (A28) 

(k — K c ) 2 



whereas the coefficient of \H a — H b \ is equal to 

2k 2 



(k - K c ) 2 



p(H - k) . (A29) 



The 2-replica correlation function G d (H a , H a ; H b , H a ) = d 2 W 2 {H a , H a ; H b , H b )/dH a dH b consists of a regular 
part and a singular part proportional to S(H a — H b ), 

6 d {H a , H a - H b , H b ) = 6 r d e9 (H a , H a - H b , H b ) + 6 s ™ 9 {H a , H a ;H a , H b )S(H a - H b ) . (A30) 

We find 

G r d e9 (H a ,H a ; H b , H b ) + m(H a , H a )m{H b , H b ) = \l + k[p(H a + k) + p(H a - k) + p(H b + k) + p(H b - k)] 

[k - K c ] 2 L 

' T " In (2 cosh -^r-) [p(H b + k) - p(H b - k)} 



k-K r k~K, 
H b , - , kH a 



In (2 cosh — ) [p(H a + k) - p(H a - k)] (A31) 



k - K c v k — K, 



2G 



and 



G s d m3 (H a , H a ;H a , H b ) = [^=V + In (2 cosh -^-)] [-^- + In (2 cosh J^-)]p{H a + k) 



k-K, 
kH a 



k-K r k-K r 



k - K r 



kff a kH b kH b 

ln(2cosh- ^)][- ^-ln(2cosh- )] p (H a - k) . (A32) 

k — K c k — K c k — K c k — K c 

(There are also singular contributions proportional to S(H b — H a ±2k) which can be discarded as we are only interested 
in the vicinity of H a = H b .) For H a , H b ->■ H and H a , H b ->■ -oo, we then find 



and 



G? B (H, H a ;H, H b ) + rh(H, H a )m(H, H b ) ~ [1 + 4kp(H - k)} 



Ah.2 fra frb 

G s ; n9 (H,H a ;H,H b ) ~ p(H - k). 

d (k-K c ) 2 ■ 



(A33) 



(A34) 



Finally, the 2-replica correlation function G d (H a ,H a ; H b , H a ) = d 2 W 2 (H a , H a ; H b , H b )/dH a dH b is given by 
G d (H a , H a ;H b , H b ) + m(H a , H a )m(H b , H b ) = 



kH a 



(k - K c ) 2 
kp(H a + k) 

+ k-Kr. 



H b — k kH b kfj b 

+ [1 + tanh ^}T(H b + k) + [1 - tanh ^}V(H b - k) 

k k-K r k-K r 



kH a , , , kH a 
In (2 cosh 



k-Kr 



kH b 

[1 - tanh ^}9(H b - H a - 2k) 

k - K c 



k-K c 

kp(H a - k) 



k k-K P 



k-K c 



kH a , , , kH a 
— In 2 cosh ■ 



H b -H a 



hfjb I. frb 

[1 - tanh ^}9(H b - H a ) + [1 + tanh ^}9(H b - H a + 2k) 

k — K r k — K r 



(A35) 



This function has a step-discontinuity at H a = H b = H (and also at H a — H b ± 2k) and therefore one must fix 
the value of 9(0) to lift the ambiguity when the sources are equal. This is done by imposing the exact symmetry 
G d (-H, -H; -H, -H) = G d (H, H; H, H) which yields 9(0) = 1/2. As a result, one has 



so that 



G d (H, H; H, H) + m(H, H)m(H, H) 



kH 



(k - K c ) 2 

kp(H + k) kH 



H — k kH kH 

— h [1 + tanh ^}V(H + k) + [1 - tanh —]V(H - k) 

k k-Kr k-Kr 



2(k-K c ) k- K c 
kp(H-k) kH 



. , kH Nr , kH 

In (2 cosh — ) | — 3 + tanh 



2(k-K c ) k-K c 



In (2 cosh 



k-Kc 
kH , 
k-K c 



} [3 + tanh 



k-K c 
kH 

k-K/ 



G d (H, H; H,H)~H 2kKH ^ [2k - H + 2kV(H - k)\ 

(k - K c y 



(A36) 



(A37) 



for H — > — oo. 

Finally, we consider the two derivatives that are needed in Eq. (96) to compute dG d (H a ; H b ) /d(H a — H b )\ Q+ in 
the single-site effective model in the limit H a , H b — > —oo. The first one is simply given by Eq. (A29), 



dG d 
d(H a - H b ) 



K c ,K d ;0+ 



2k 2 
[k-K c 



-p(H - k) 



(A38) 
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The second one is obtained from Eq. ( A33 1 , 



dG d 



lim 



— [G d (H a ,H a ; H\ H b ) + m(H a , H a )fn(H b , H b )} 



dK d {H a -H b ) o+ 



— OO 



dH a dH b 



1 + 4kp(H - k) 
(k - K c ) 2 



(A39) 



where we recall that the subscript + indicates the limit of equal fields, H a = H b = H . 
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